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An  approach  to  multivariate  interpolation  is  described.  The  algorithm  is 
applicable  in  arbitrary  dimensions,  and  can  generate  surfaces  of  arbitrary 
smoothness.  This  is  accomplished  by  tesselating  the  (polyhedral)  domain  into 
simplices  and  using  one  dimensional  algorithms  to  construct  interpolants  first 
on  edges  and  then  successively  on  higher  order  faces  by  blending  methods.  The 
result  is  a  piecewise  rational  function  of  a  high  degree  which  has  the 
prescribed  global  smoothness  and  interpolates  to  the  original  data.  The 
interpolants  are  local,  i.e.  their  evaluation  at  a  point  requires  only  data  on 
the  simplex  that  the  point  resides  in.  The  schemes  require  data  of  the  same 
degree  as  the  degree  of  global  smoothness.  The  degree  of  polynomial  precision 
is  greater  than  or  equal  to  the  degree  of  smoothness.  The  approach  derives 
its  power  and  simplicity  from  the  fact  that  derivatives  in  directions 
perpendicularly  across  faces  are  incorporated  directly  as  data.  .. 
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SIGNIFICANCE  AND  EXPLANATION 


An  approach  is  described  to  the  interpolation  of  data  in  arbitrarily  many 
variables.  The  schemes  are  of  arbitrary  smoothness  and  require  data  of  the 
same  degree  as  the  degree  of  smoothness.  The  domain  is  assumed  to  be 
tesselated  into  simplices  (e.g.  triangulated  in  the  case  of  two  variables). 
Evaluation  at  a  point  requires  only  data  on  the  simplex  that  the  point  resides 
in.  The  methods  described  here  constitute  the  only  known  local  interpolation 
schemes  in  arbitrarily  many  variables  and  with  an  arbitrary  degree  of  smooth¬ 
ness. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


MULTIVARIATE  PERPENDICULAR  INTERPOLATION 


Peter  Alfeld* 

1 .  Intr oduc t ion 

The  problem  of  interpolating  to  scattered  multivariate  data  is  becoming  increasingly 
important.  Applications  include  the  modeling  of  physical  phenomena  involving  space  and 
time  (e.g.  combustion,  temperature,  pressure,  etc.)  and  the  design  of  geometric  objects 
(e.g.  the  body  of  a  car  or  an  aircraft).  For  a  recent  survey  of  this  area  see  Barnhill, 
1983. 

In  this  paper,  we  describe  an  approach  to  multivariate  interpolation  in  an  arbitrary 
number  of  variables,  and  with  an  arbitrary  degree  of  smoothness.  He  assume  that  the  domain 
of  interest  has  been  tesselated  into  simplices,  and  that  at  the  data  points  we  are  given 
the  values  of  all  partial  derivatives  through  the  same  order  as  the  degree  of  smoothness. 
Any  other  directional  derivatives  can  then  be  constructed  from  the  partial  derivatives.  We 
do  not  address  the  question  of  how  the  tesselation  may  be  accomplished,  or  how  higher  order 
data  may  be  generated  from  lower  order  data.  For  some  answers  to  these  questions  see 
Barnhill  and  Little,  1983,  or  Alfeld,  1984a. 

Our  interpolants  are  piecewise  rational  functions  of  the  prescribed  smoothness  that 
Interpolate  to  the  given  data.  The  basic  idea  consists  of  constructing  the  interpolant  on 
an  individual  simplex  as  a  blend  of  lower  dimensional  interpolants  on  faces  of  the 
simplex.  Recursion  of  this  procedure  ultimately  leads  to  one  dimensional  problems  which 
are  easy  to  solve. 

The  interpolants  derive  their  power  and  simplicity  from  the  fact  that  derivatives  in 
directions  perpendicular  to  faces  of  a  simplex  are  incorporated  directly  as  data.  Such 
derivatives  govern  the  smoothness  of  the  interpolants  between  simplices.  If  derivatives  in 
other  directions,  e.g.  parallel  to  edges,  are  used,  then  an  algebraic  link  to  perpendicular 
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derivatives  has  to  be  constructed.  This  leads  to  extremely  complex  algebraic  manipulations 
that,  moreover,  at  the  present  state  of  the  art,  cannot  be  carried  out  for  a  general  degree 
of  smoothness  and  a  general  number  of  variables.  The  central  role  of  perpendicular 
derivatives  gave  rise  to  the  title  of  this  paper. 

The  schemes  proposed  here  are  local,  i.e.  their  evaluation  at  a  point  requires  only 
data  on  the  simplex  that  the  point  resides  in.  In  two  variables  there  are  many  local  C1 
and  c  schemes  (see  Barnhill,  1983,  and  the  references  therein).  There  are  very  few 
local  trivariate  C1  schemes  (see  Alfeld,  1984,  and  Barnhill  and  Little,  1983)  and  no 
alternative  c  schemes.  The  methods  proposed  here  are  the  only  existing  local  schemes 
for  scattered  data  in  an  arbitrary  number  of  variables  and  with  an  arbitrary  degree  of 
smoothness. 

The  paper  is  organized  as  follows:  In  section  2,  we  develop  some  geometric  concepts 
and  notation.  In  section  3,  we  described  two  types  of  perpendicular  interpolants.  In 
section  4,  some  formulas  are  given  that  are  useful  for  the  implementation  of  bi-  and  tri¬ 
variate  schemes;  and  section  5  contains  (trivariate)  numerical  examples  and  comparisons. 

2.  The  Geometry  of  Perpendicular  Interpolants 

We  shall  he  ultimately  concerned  with  polyhedral  regions  in  tP  that  have  been 

tesselated  into  n-dimensional  simplices.  Any  such  simplex  is  defined  by  its  vertices  which 

we  denote  by  v1#  v2'"’,vn+1‘  *  general  point  x  e  rf1  is  described  in  terms  of 

barycentrlc  coordinates  b1,...,bn+1  as 

n+ 1  n+ 1 

x  -  l  b.V,  where  l  b,  =  1  .  (2.1) 

i=1  i-1 

We  will  always  assume  that  the  general  simplex  is  non-degenerate,  i.e.  that  the  linear 
system  (2.1)  possesses  a  unique  solution. 

Most  of  the  sequel  will  be  concerned  with  the  development  of  the  interpolants  as  they 
are  defined  in  the  context  of  a  single  simplex.  It  will  be  convenient  to  think  of  the 
simplex  as  the  affine  space  spanned  by  the  vertices.  We  define 


(2.3) 


bound(S)  (x  e  S  I  3  1  e  { 1,2, . . . ,n+l}  s  *  Ot  . 

When  we  need  to  refer  to  the  simplex  as  the  convex  hull  of  its  vertices  we  will  use 
the  notation 

n+1  n 

conv(S)  »  {x  =  }  b  V  |  \  b.  »  1,  b  >  0  i  -  1,...,n+l} 

i=1  i=1 

We  will  denote  by  F  that  n-1  dimensional  face  of  S  that  is  obtained  by  removing 
the  vertex  V^.  F  is  itself  a  simplex  and  the  notations  bound  (F)  and  conv  (F)  apply. 

The  set  of  all  facets  of  S  of  dimension  p,  say,  is  denoted  by  Fy/  and  the  set  of  all 
facets  is  F. 

The  anchor  of  F  is  that  face  of  S  whose  normal  forms  the  smallest  angle  with  the 
normal  of  F.  With  the  point  x  and  the  face  F  we  associate  the  straight  line  through 
x  perpendicular  to  F.  That  line  is  the  line  of  fixation  of  x  with  respect  to  F.  F 
is  also  called  the  base  face  of  the  line  of  fixation.  The  point  B^tx)  where  the  line  of 
fixation  intersects  its  base  face  is  the  base  point  of  x  and  the  point  Ap(x)  where  it 
intersects  the  anchor  is  the  top  point  of  x  (all  with  respect  to  the  base  face). 

Before  going  into  algebraic  details  in  section  3  it  is  useful  to  visualise  perpendicu¬ 
lar  interpolants  geometrically.  Consider  a  fixed  point  x  and  a  fixed  simplex  S.  In 
unanchored  interpolants,  we  interpolate  along  all  lines  of  fixation  through  x  to  all 
derivatives  through  given  order  at  the  base  points.  The  final  value  of  the  interpolant  is 
obtained  by  suitably  blending  the  values  of  the  univariate  interpolants  defined  on  the 
lines  of  fixation.  In  anchored  interpolants,  we  interpolate  additionally  to  function  (but 
no  derivative)  values  at  the  top  points  of  x. 

Note  that  perpendicular  interpolants  require  data  on  the  faces  of  S.  These  are  in 
turn  defined  in  terms  of  lower  dimensional  interpolants,  giving  rise  to  the  recursive 
structure  of  perpendicular  interpolants.  It  is  possible  that  the  top  point  of  the  line  of 
fixation  lies  outside  the  convex  hull  of  the  points  defining  the  anchor,  necessitating 
extrapolation.  This  is  why  in  definition  (2.2)  and  (2.3)  we  do  not  restrict  ourselves  to 
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the  convex  hull  of  {v^}.  In  t^ie  presence  of  obtuse  angles  extrapolation  may  also  be 
necessary  in  the  base  face. 

The  above  geometric  concepts  are  illustrated  in  Figures  1  and  2  for  the  case  that 
n  «  2.  In  section  4  on  computational  aspects  algebraic  expressions  are  given  for  n  =  2 
and  n  •  3. 


3.  The  Interpolants 

3.1.  Uhanchored  Interpolants 

We  shall  construct  operators  P^  which  associate  to  each  f  a  function  Pgf  S  Cm 
which  interpolates  f  and  its  derivatives  of  order  <  m  at  the  vertices  of  S. 

Actually  P^f  only  depends  on  the  values  of  f  and  its  derivative  values  at  the  vertices 
and  hence  in  actuality  is  defined  for  each  set  of  discrete  data.  However,  it  is  convenient 
to  think  of  P^  as  operating  on  functions.  If  P^  is  applied  piecewise  on  a  tesselation 
of  a  polyhedral  domain  D  into  simplices,  then  the  resulting  piecewise  defined  function  is 
in  (^"(D)  • 

Now  the  operator  P^  is  built  up  from  the  lower  dimensional  operators  p£  where  F 
is  any  face  of  S  of  dimension  n  -  1  and  0  <  k-  <  m.  In  this  way,  the  operators  P^J 
are  defined  recursively  starting  with  simplices  of  dimension  0  (the  vertices)  and  working 


up  to  higher  dimensional  simplices. 

m  |r 

The  operator  Pg  is  gotten  by  blending  the  lower  dimensional  operators  Pp. 

Suppose  is  a  vertex  of  S  and  F  is  the  n  -  1  dimensional  face  of  S  which  does 

not  contain  V^.  For  any  point  x  as  in  (2.1)  which  is  not  in  bound  (S)  we  define 


The  exponents  of  the  barycentric  coordinates  br  are  chosen  to  be  even  so  that  the 


blending  functions  Cp  constitute  a  nonnegative  partition  of  unity  even  if  some  bary- 
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Figure  1:  Extrapolation  in  the  anchor 
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centric  coordinates  are  negative.  This  is  necessary  because  as  pointed  out  in  section  2 
some  extrapolation  on  faces  of  S  may  be  necessary. 


We  do  not  define  c®  for  points  in  bound(S)j  however,  since  b^  =  0  on  F,  c®  is  to 

be  thought  of  as  a  function  which  is  1  on  F  and  0  on  the  other  faces  of  S.  For 

such  a  face  F  we  also  introduce  the  operators  Gp  defined  by 

v 


G™(f)Cx)  :=  “  ^PF_V(1~f)  (VX)) 
v=0  asF 


(3.2) 


where  the  p£  are  the  lower  dimensional  operators,  and,  as  introduced  earlier,  Bp(x)  is 
the  base  point  of  x  and  sF  is  the  normal  to  F,  normalized  to  have  length  equal  to  the 
distance  from  to  F.  Note  that  the  operator  G®  extends  data  on  F  by  Taylor 

interpolation  in  the  direction  perpendicular  to  F,  hence  the  name  perpendicular 
interpolation . 

Then,  for  n  ■*  dim(S)  >  1,  and  m  >  0,  the  operators  Pg  are  defined  by 


p”(f)(x) 


(  1  c“(x)G”(f ) (x) 

I  FeFn 


x  t  bound ( S ) 
x  e  f  e  F  , 

n-1 


(3.3) 


We  shall  see  shortly  that  P^(f)(x)  is  well  defined  on  bound(S).  Before  doing  this, 
however,  we  must  define  P^  if  dim(S)  <1  or  m  =  0.  This  is  accomplished  by 


n+1  n+1 

P°<n  (  l  b  v  )  :=  )  b  f(V  )  (3.4) 

S  i=1  1  1  i=1  1  1 

and 


dim(S)  =  0  =>  P™(f)(x)  :=  f(V)  (S  -  {v}  )  (3.5) 

s 

dim(S)  =  1  ==> 


p”(f)(tv2  +  (1-t)Vt) 


m 

-i=n 


rv,™  in  jig 
'  0,j(t)  3(V  -V.) 


<V 


hm  ...  ajf 

h1rj  3(V,-V.) 


(v2)i 


(3.6) 


where  S  is  spanned  hy  V1  and  Vj  and  the  h®  ^  are  cardinal  polynomials  of  degree 
2m  +  1  that  satisfy  and  are  defined  uniquely  by  the  properties 

=  5  S.. 

su  xk 

t=k 

s,li  -  0,1,.. .,m;  i,k  -  0,1,  and  5  the  Kronecker  delta. 

Thus  (3.4)  describes  multivariate  linear  interpolation  and  (3.6)  univariate  Hermite 
Interpolation  of  degree  2m  +  1. 

We  turn  now  to  the  smoothness  of  P^f •  ’  s  clear  that  P^f  is  infinitely  often 

differentiable  at  all  points  of  S  which  a  '<it  in  bound (S).  To  study  P®f  on  bound(S) 

we  need  the  following  lemma: 

Lemma.  The  blending  functions  c®  have  the  following  properties: 

1.  For  any  differentiation  operator  p  of  order  t,  1  <  l  <  m,  and  any  Xg  e  bound(S) 

lim  pc®(x)  *  0 

^Xo 

xflbound  ( S ) 

2.  With  G  »  (F  |  xQ  e  F)  /  g 

lim  )  c”(x)  =  1 

x*Xg  FeG 

x^bound(s) 

Proof:  The  lenna  is  proved  by  inspection. 

We  now  come  to  the  fundamental  theorem  of  this  section: 

Theorem  1 : 

The  operator  P^  has  the  properties: 

(i)  P^f  is  well  defined  for  each  f. 

(ii)  p£f  e  C®  for  each  f. 
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(iii)  For  any  Xq  e  F  e  Fn_-j  and  any  differential  operator  V  of  order  I  <  m  involving 
only  directions  contained  in  S: 

PPg<fXx0>  =  P^k(Pf)(x0)  . 


Proof:  The  theorem  is  proved  by  induction  on  dim(S).  The  statements  are  clear  from  the 

definition  (3.6)  of  P^  if  dim  (s)  =  1.  Now  suppose  we  have  established  the  theorem  for 

all  simplices  of  dimension  n  -  1  and  let  S  be  a  simplex  of  dimension  n. 

Consider  first  statement  (i).  If  xQ  e  bound(s)  and  xQ  e  F  n  F’  for  two  faces  F 

and  F'  of  S  then  Xq  e  G  with  G  a  facet  of  dimension  n  -  2  contained  in  both  F 

and  F' .  But  then  P^!(f)(x0)  :=  PG(f)(Xg)  and  similarly  for  F'  by  the  induction 
hypothesis.  This  establishes  (i). 

Now  consider  (iii)  and  suppose  F  e  Fn-1  a  fixed  face  of  S.  Without  loss  of 
generality  we  may  assume  that  V  has  been  written  as  the  sum  of  operators  of  the  form 
OpD  where  is  of  order  j  involving  only  derivatives  in  directions  contained  in 

3kf 

F  and  V  =  — —  with  S_  the  normal  to  F  and  j+k  =  l  <  m.  We  first  observe  that 
3sk 

F 

lim  P(GFf)(x>  -  p”-*(Pf)(x0>  .  (3.7) 

*+x0 

x^F 
—  V 

Indeed,  Ppb^  *  0/  and  P(b^)  =0  at  unless  v  =  k  (in  which  case 

—  k 

Pfb^)  *  k!).  Hence,  from  our  induction  hypothesis 


m  __  D,  v 

lim  V(G  f )  ( x)  =  lim  l  (Pe))(PpPrVt^MB  (x,)) 


„vvJ'kF  F  \  V  F 

**xo  v=i  asF 


lim  t»F(p”"k(^-|)](BF(x)' 
**xo  38f 


-  pf'*(Pf  7T><VV> 

SSF 


=  P™’4(Pf)(x0)  . 


-9- 


This  shows  (3.7) 


The  letter  equality  holds  because  P 


akf 

F  3sk 

F 


Pf  and  Bp(  Xg )  «  xQ. 


From  (3.7)  we  can  conclude  that  under  the  same  hypotheses 

lim  P(Pgf )  (x)  -  p“‘*(0f)(xo)  . 

x+x0 


(3.8) 


x$F 

To  see  this,  we  first  observe  that  as  in  the  proof  of  (i),  the  right  hand  side  in  (3.8)  is 
well  defined  if  Xq  lies  in  two  different  n  -  1  dimensional  faces  of  S.  Let  us  denote 
by  L  the  common  value  of  pST^PfHXg)  for  all  faces  F'  that  contain  xQ.  Using  the 


Lemma  and  (3.7)  we  have 


lim  P(P™f)(x)  =  lim  (  c™(x))l  «  L 


(3.9) 


x*x„ 


x^bound ( S ) 


x»x0  F 

x0eF 


The  last  equality  holds  because  in  the  limit  all  derivatives  of  all  c”,  and  all 

c™  with  Xg  ^  F  vanish,  rendering  the  sum  in  (3.9)  equal  to  1,  and  because  the  Gpf(x) 

are  continuous  and  hence  bounded  in  a  neighborhood  of  xQ.  This  establishes  (iii). 

The  statement  (ii)  now  follows  by  induction  on  the  degree  £  <  m  of  a  differentiation 

m 

operator.  Let  x0  e  bound(S).  In  the  proof  of  (iii)  we  have  seen  that  lira  Pgf(x)  is 

x*x0 

well-defined  and  assumes  the  value  required  in  (3.3).  It  is  hence  continuous.  Now  assume 
1-1 

P^  e  C  (S),  £  <  m,  and  P  is  a  differentiation  operator  of  degree  £ .  Then 

lim  Pp”f(x)  is  again  well  defined.  Since  P™  8  C*  (S),  Pp*J  is  continuous.  This 
s  s  s 

x*x0 

establishes  (ii)  and  completes  the  proof  of  the  theorem. 


The  properties  of  interpolation  and  global  smoothness  of  piecewise  applied  unanchored 
interpolants  follow  immediately  from  Theorem  1. 


Corollary:  Assume  a  polyhedral  region  D  has  been  tesselated  into  simplices  and  let 
P®(f)(x)  =  Pg(f)(x)  whenever  x  e  conv(S).  Then  Pmf  is  well  defined,  P-mf  e  (^(D), 
and  for  all  vertices  V  and  differential  operators  P  of  order  k  <  m 

P  (Pmf)(V)  =  Pf(V)  . 
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Proof:  Nothing  needs  to  be  shown  in  the  interior  of  a  simplex.  So  let  P  be  a 
comnon  idcet  of  two  simplices  S  and  S'.  Then  by  a  possibly  repeated  application  of 
Theorem  1  (iii). 


Hence  Pmf  and  its  derivatives  are  well  defined.  Its  smoothness  follows  as  in  the  proof 
of  theorem  1  (ii).  Interpolation  follows  from  (3.10)  and  the  definition  (3.5). 

We  now  turn  to  the  question  of  polynomial  precision. 

Theorem  2:  Let  f  be  any  polynomial  in  n  variables  of  degree  <  m,  and  let  P®f 

be  defined  as  in  the  corollary  of  Theorem  1.  Then 

P^f  -  f  . 

Proof:  It  is  sufficient  to  consider  only  an  individual  simplex  S.  If  m  «  0,  the 

operator  P^  even  reproduces  linear  functions  which  are  of  degree  m  +  1.  For  m  >  0, 

the  proof  is  again  by  induction  on  dim  (S).  The  statement  is  clear  from  the  definition 
(3.6)  if  dim  (s)  -  1  since  univariate  Hermite  Interpolation  reproduces  polynomials  of 
degree  up  to  2m  +  1.  This  follows  from  the  existence  and  uniqueness  of  the  interpolant 
which  is  shown  e.g.  in  Davis,  1975,  p.  29. 

Now  suppose  we  have  established  the  theorem  for  all  simplices  of  dimension  n  -  1  and 
let  S  be  a  simplex  of  dimension  n.  Let  f  be  a  polynomial  of  degree  <  m.  Then  f 
reduces  to  a  univariate  polynomial  of  degree  <  m  along  any  line  of  fixation.  It  follows 
that  for  any  face  F  of  S  the  operator  G®  reproduces  f  exactly  because  by  (3.2) 

Gp  just  applies  Taylor  interpolation  to  data  on  F  which  are  exact  by  the  induction 
hypothesis.  The  theorem  now  follows  from  (3.3)  since  the  c®  constitute  a  partition  of 


unity. 


,%V* 

\*/ 


\\‘J 


3.2  Anchored  Interpolants 

In  this  subsection,  we  increase  the  degree  of  precision  of  unanchored  interpolants  by 
also  interpolating  to  function  (but  no  derivative)  values  at  the  top  point  of  the  line  of 
fixation.  It  is  infeasible  to  interpolate  to  derivatives  at  the  top  point  since  the 
direction  of  interpolation  is  not  perpendicular  to  the  anchor.  The  necessary  changes  in 
the  interpolant  consist  of  the  addition  of  another  term  in  the  definition  of  the  face 
operators  and  a  modification  of  the  cardinal  functions.  The  notation  is  essentially  the 
same  as  that  for  unanchored  interpolants,  except  that  differing  objects  have  been 
distinguished  by  bars.  We  define 


^<x> !’  v-o  ^(bp,tp)  ^r^)(vx,) 

p 


+  WVV 


where  A  is  the  anchor  of  F,  TF(x)  is  the  top  point  of  the  line  of  fixation  with  respect 
to  F  through  x,  and  Bj.(x)  and  sF  are  as  before.  The  quantities  bp  and  tp  are 
defined  by 

bF  lx  -  Bf(x)«2 

and 


tp  lx  -  Tp(x) I 2 
— m 

and  the  pv  replace  the  Taylor  polynomials.  They  are  polynomials  of  degree  m  +  1 
defined  by  the  cardinal  properties 


3r  p^(b,t) 
3br 


b=0 


r  =  0, 1 , . . .  ,m 
v  =  0,1,... ,m+1 


and  P?<t,t)  -  «rm+1  r  ”  0,1, 
Explicitly  they  are  given  by 


,m+1 


m,v  b_  ,  b  m-v+1. 

P„(b't)  "  ~  ^  ‘  (v>  > 
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A. 

y. 


*-.v 


,v. 


V.V 


*  - '  * « •  •  ’ 

*  i  '  •  »  1 


thus  the  operator  G  extends  data  on  F  and  its  anchor  by  univariate  interpolation  along 
F 

the  line  of  fixation. 

— m  _n» 

The  operators  P_  are  now  defined  analogously  to  vz  by 
S  o 

x  ^  bound (S) 

x  e  F  e  F  . 

n-1 

if  dim  (S)  >  1  and  m  >  0,  and 

*>>  ■  >"<*>'  ■  >s(x> 

otherwise. 

The  following  theorem  is  analogous  to  Theorem  1 : 


re  F 


Cp(x  )G“(  f )  (x) 


P“  (f)(x> 
s 


n- 1 
PjifK*) 


Theorem  3i 

— ni 

The  operator  Pg  has  the  properties 

(i)  p"f  is  well  defined  for  each  f 
s 

(li)  ^f  e  c"  for  each  F 

(iii)  FOr  any  ij  !  F  6  ^n-1  an^  “"F  differential  operator  V  of  order  t  <  m 
Involving  only  directions  contained  in  Ss 

tfJ'fMV  -^*<P«*X0>  • 

Proofs  The  proof  of  Theorem  is  like  that  of  Theorem  1.  The  only  difference  is  that 

the  G^(f )  differ  from  the  G^!(t)  on  faces  other  than  F.  However,  all  that  is  required 
-m 

is  the  boundedness  of  G^IfHx)  as  x  approaches  a  point  in  a  face  other  than  F.  This, 
however,  follows  from  the  continuity  of  G^(f). 

Global  Smoothness  and  interpolation  follows  as  for  Theorem  1: 


Corollary:  Assume  a  polyhedral  region  D  has  been  tesselated  into  sioplices  and  let 

T'lfllx)  -  7"<f)(x)  whenever  x  e  oonv(S).  Then  P^f  is  well  defined,  P^f  e  Cm(D),  and 
s 

for  all  vertices  V  and  differential  operators  V  of  order  k  <  m 

0(p”f)(V)  -  0f(V)  . 

The  advantage  of  anchored  interpolants  is  that  the  degree  of  polynomial  precision  is 
increased  without  raising  the  degree  of  the  required  data: 

Theorem  4:  Let  f  he  any  polynomial  in  n  variables  of  degree  <  m+1,  and  let 
“in 

P  f  be  defined  as  in  the  corollary  of  Theorem  3.  Then 

i^f  -  f  . 

The  proof  of  theorem  4  is  analogous  to  that  of  theorem  2. 


4.  Implementation  of  Interpolants  in  Two  or  Three  Variables 

In  this  section,  we  develop  algebraic  representations  of  the  geometric  concepts 
introduced  in  the  preceding  sections.  We  restrict  our  attention  to  two  or  three 


variables.  An  edge  of  a  simplex  will  be  denoted  by  e 


ij 


V 


4.1  Two  Variables 

In  this  case,  the  generic  simplex  S  is  a  triangle,  and  its  faces  are  edges.  The 
anchor  A  of  the  base  face  e  is  that  edge  other  than  e  that  maximizes  the  expression 


| cos  a | 


I  eTA  I 


Iel2 IAIj 


where  a  is  the  angle  formed  by  e  and  A.  In  the  case  of  an  isosceles  triangle,  the 
resulting  tie  may  be  broken  arbitrarily,  but  this  must  be  done  consistently.  Some  care  may 
be  necessary  to  avoid  different  assignments  of  the  anchor  at  different  stages  of  the 
computation  due  to  round-off  errors. 


To  maintain  generality  we  consider  a  general  triangle  S  with  vertices  Vj,  V^,  Vk 


where  e^  is  the  base  face  and  e^j  is  its  anchor,  as  depicted  in  Figures  1  and  2.  Any 
actually  occuring  situation  can  be  described  by  assigning  suitable  values  to  i,  j,  and 
k. 

The  normal  4  to  the  base  face  is  expressed  in  the  form 


e,  .  +  ye  .  where 

1  j  jK 


.T 

♦  e 


jk 


0  i.e. 


T 

*ij*jk 


(4.1) 


T 

®jkejk 


A  general  point  x  is  written  in  terras  of  barycentric  coordinates  as 
x  -  biVi  +  bjVj  +  Vk'  bi  +  bj  +  bk  *  1  • 


Let  B  denote  the  base  point  of  the  line  of  fixation  through  x,  and  let  T  denote 


its  top  point,  with  the  normalization  implied  by  (4.1)  we  have  that  h,  is  the  natural 

3bi 

parameter  along  ♦,  i.e.  ”  1.  we  obtain 


x  -  B  t  b,*,  B  -  ajVj  +  akVk 
where  aj  ”  bj  +  “  1  ~  and 


T  ”  civi.  +  cjvj  where 


—  and  c,  *  1-c, 
y  j  i 


Note  that  always  y  ?  0  because  the  choice  of  the  anchor  avoids  right  angles  o. 


4.2.  Three  Variables 

We  now  repeat  this  development  for  the  case  of  three  variables  and  consider  a  general 
tetrahedron  S  with  vertices  Vi#  Vj,  Vk,  Vj  where  the  base  face  is  spanned  by  Vj,  Vk, 


and  Vj.  The  normal  to  the  base  face  is  written  as 


♦  “  eij  +  aeJk  *  Pek* 


where  a  and  p  are  defined  by  the  linear  system 


*Te  >  *Te  -  0  . 

jk  kt 
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Explicit  expressions  for  a  ana  8  are  given  by: 


T  T  T  1 

*ikekt  i1®kt  "  ktek*ei 


*jkejkektekt  ~  6 jk®kt ' 


kekiei1®1k  eikeikeilekt 


ejkejkektekt  "  <ejkekt) 


The  anchor  of  the  base  face  is  that  face  other  than  the  base  face  whose  normal  q 


say,  maximizes  the  expression 


,®,2,q,2  ’ 


A  general  point  x  is  now  written  as 


X  -  biVi  ♦  bjVj  ♦  bkvk  ♦  btvt 


bi +  bj +  bk +  bl  • 1  • 


Again  we  denote  the  base  point  by  B  and  the  top  point  by  T.  We  obtain 


8  *  Vi  *  Vk  *  Vi 


“j  "  bj  +  (1-°,bi 


■  *V  +  (o-B  )b. 


*k  *  bk 


at  *  1  “  aj  “  ak 


and,  assuming  the  anchor  is  spanned  by  V^,  V ^ ,  Vk 


T  -  c^  +  ciVi  *  ckVk 


'i  8 


Cj  *  a^  +  (n-DCj^ 


ck  “  1  -  ci  “  cj  ‘ 


5.  Computational  Aspects  and  Numerical  Results 
5.1.  Computational  Aspects 

He  have  Implemented  all  the  schemes  described  In  this  paper  Into  FORTRAN  research 
codes.  A  major  tool  In  constructing  the  codes  was  the  symbol  manipulation  language  RSDOCB 
(Hearn,  1983),  which  can  be  instructed  to  write  mathematical  formulas  directly  in  PomiUUi 
notation.  He  have  examined  all  interpolants  using  the  smoothness  tester  MXCnOBCOHt  (Alfeld 
and  Harris,  1984).  Hiis  investigation  confirmed  that  the  codes  do  indeed  possess  the 
smoothness  and  precision  properties  implied  by  the  mathematical  construction. 


5.2.  Numerical  Results 

In  this  section,  we  compare  several  trivariate  schemes  on  a  simple  problem.  Hie 
domain  is  a  distorted  cube  that  has  been  tesselated  into  12  tetrahedra.  The  tesselation  is 
described  in  detail  in  Alfeld,  1984b.  Hie  guiding  principle  behind  the  construction  of  the 
domain  was  to  avoid  artifacts  due  to  e.g.  edges  or  faces  being  parallel  to  coordinate  axes. 


For  ease  of  comparison  an  underlying  function  was  used.  This  was  given  by 

F(x,y,s)  «  /r  -X  -y  -z  (5.1 

for  several  values  of  r.  Hie  derivative  data  required  for  the  schemes  were  exact.  Hie 
distance  of  the  vertices  of  the  domain  from  the  origin  varied  between  18.71  and  187.6. 

For  comparison,  we  also  include  transfinite  variants  of  our  interpolants,  i.e., 
schemes  where  the  data  on  faces  are  given  exactly  rather  than  being  defined  by  lower 

■  n 

dimensional  interpolants.  Thus  w  and  w  denote  the  transfinite  interpolants 
corresponding  to  P*"  and  P® ,  respectively. 

Table  1  contains  the  numerical  results.  Hie  interpolation  schemes  are  ordered  by 
increasing  complexity  as  defined  by  the  amount  of  CPU  time  (on  the  DEC  System  20  of  the 
College  of  Science  at  the  University  of  Utah)  required  for  one  evaluation  of  the 


interpolant.  The  individual  schemes  were  not  programmed  in  the  most  efficient  fashion 


possible,  but  their  implementations  are  sufficiently  similar  such  that  the  comparison  is 
valid.  The  following  is  a  more  detailed  description  of  the  table  by  columns: 

1.  (Scheme)  Short  notation  for  the  scheme.  For  most  of  the  scheme  the  notation  is 
defined  in  this  paper,  but  the  schemes  defined  in  Alfeld,  1984b  are  included  for 
comparison. 

2.  (Discrete  versus  Transfinite)  Indicates  whether  the  scheme  is  a  transfinite  or  a 
discrete  one. 

3.  (Anchor)  Indicates  whether  the  scheme  is  anchored  or  not.  N.  A.  means  that  the 
scheme  is  not  based  on  perpendicular  interpolation  so  that  this  distinction  does  not  apply. 

4.  (Precision)  All  (trivariate)  polynomial  of  degree  up  to  the  given  number  are 
reproduced  exactly  by  the  scheme. 

5.  (Smoothness)  All  derivatives  through  the  given  order  are  continuous. 

6.  (r)  The  value  of  r  in  (5.1). 

7.  (Error)  The  maximum  relative  error  (computed  over  the  entire  distorted  cube). 

This  quantity  was  computed  by  searching  in  each  tetrahedron  with  Winfields  Method 
(Winfield,  1973)  starting  at  the  centroid  of  the  tetrahedrom.  This  algorithm  identifies 
local  maxima.  Hence  there  is  a  slim  possibility  that  the  overall  global  maximum  was 
overlooked  by  the  method.  However,  this  is  an  unlikely  occurrence  because  on  each  tetra¬ 
hedron  the  primitive  function  (although  not  the  relative  error)  is  uniformly  convex. 

8.  (CPU)  The  CPU  time  in  milliseconds  for  one  evaluation  of  the  interpolant  (not 
including  the  setting  up  of  the  relevant  data  structure). 

9.  (COMM.)  Comments: 

(1)  This  is  the  standard  linear  interpolant  that  is  included  here  for  comparison. 

(2)  This  is  the  transfinite  scheme  described  in  Alfeld,  1984b.  It  is  also  included 
for  comparison. 

(3)  Here  the  necessary  extrapolation  on  the  anchor  caused  an  evaluation  of  the 

2  2  2  2 

primitive  function  outside  of  its  domain  where  x  +  y  +  z  >  r  .  The  computing  system 
proceeded  by  replacing  negative  radicands  by  their  absolute  values. 
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(4)  Thia  is  the  discrete  interpolant  described  in  Alfeld,  1984b.  It  is  also  included 
for  comparison. 

The  following  points  emerge  from  the  table: 

1.  Most  schemes  give  reasonable  results ,  considering  that  a  large  part  of  the  domain  of 
the  primitive  function  is  covered  by  only  12  tetrahedra. 

2.  The  discrete  schemes  are  significantly  more  complex  than  the  transfinite  schemes.  The 
only  exception  to  this  is  the  very  simple  linear  interpolation  scheme,  which  yields  a 
surface  that  is  only  continuous. 

3.  The  discrete  anchored  C1  interpolant  P1  is  much  more  efficient  than  its  competitor 
described  in  Alfeld,  1984b  which  has  very  similar  properties.  Indeed,  an  attempt  was  made 
to  construct  a  C  scheme  along  the  lines  described  in  that  paper.  This  attempt  failed 

2 

for  purely  practical  reasons,  the  C  scheme  was  simply  too  complex  for  all  practical 
purposes.  Since  the  smoothness  and  precision  of  the  C1  schemes  are  identical,  and  the 
maximum  relative  errors  are  very  similar,  this  demonstrates  the  superiority  of  the 
perpendicular  interpolation  approach  described  in  this  paper. 

4.  As  one  would  expect,  the  discretization  of  a  transfinite  scheme  increases  the  error 
substantially.  This  is  a  price  that  simply  has  to  be  paid  in  a  practical  application  since 
transfinite  data  will  usually  not  be  available. 

5.  The  accuracy  also  deteriorates  as  the  radius  of  the  domain  decreases  and  the  boundary 

of  the  domain  approaches  the  outmost  point  of  the  tesselation. 

2  1 

6.  The  gain  in  accuracy  due  to  employing  a  C  rather  than  a  C  scheme  is  surprisingly 
small. 
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